Single cell RNA-seq of PBMCs were collected at baseline and stimulated with plaugue. The effect of stimulation was estimated in each celltype using pseudobulk expression. The resutling effect sizes across cell types were restimated using MASH. Differentially expressed genes in a given cell type are identified as those with \(\text{lfsr} < 0.01\).
library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
data <-new.env()load('data/plague.RData', env=data)make_list_and_background <-function(gene, lfsr){ gene_list <- gene[which(lfsr <0.01)]return(list(gene_list=gene_list, gene_background=gene))}library(AnnotationDbi)
Loading required package: stats4
Loading required package: BiocGenerics
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:dplyr':
combine, intersect, setdiff, union
The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabs
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: IRanges
Loading required package: S4Vectors
Attaching package: 'S4Vectors'
The following objects are masked from 'package:dplyr':
first, rename
The following object is masked from 'package:utils':
findMatches
The following objects are masked from 'package:base':
expand.grid, I, unname
Attaching package: 'IRanges'
The following objects are masked from 'package:dplyr':
collapse, desc, slice
Attaching package: 'AnnotationDbi'
The following object is masked from 'package:dplyr':
select
One thing to note is that the gene lists are quite large when using an lfsr threshold of 0.01, and the DE genes in each cell type have substantial overlap.
# Sample list of listslist_of_lists <- purrr::map(gene_lists$gene_list, ~.$gene_list)names(list_of_lists) <- gene_lists$celltype# Function to calculate the size of the intersection of two listsintersection_size <-function(a, b) {length(intersect(a, b))}# Determine the number of listsn <-length(list_of_lists)# Create an empty matrix to store the intersection sizesintersection_matrix <-matrix(0, n, n, dimnames =list(names(list_of_lists), names(list_of_lists)))# Calculate the size of the intersection for each pairfor (i inseq_along(list_of_lists)) {for (j inseq_along(list_of_lists)) {# Fill in the matrix with intersection sizes, ensuring i <= j to avoid redundant calculationsif (i <= j) { intersection_matrix[i, j] <-intersection_size(list_of_lists[[i]], list_of_lists[[j]]) } }}# Since the intersection operation is symmetric, mirror the upper triangle to the lower triangleintersection_matrix[lower.tri(intersection_matrix)] <-t(intersection_matrix)[lower.tri(intersection_matrix)]# Display the matrixpheatmap::pheatmap(intersection_matrix, cluster_rows =FALSE, cluster_cols =FALSE, display_numbers =TRUE)
# Function to find unique elements in a list compared to othersfind_unique_elements <-function(target_list, other_lists) {# Combine all other lists into a single vector combined_others <-unlist(other_lists)# Find elements in the target list that are not in the combined vector of other listssetdiff(target_list, combined_others)}# Create a new list to store lists of unique elementsunique_elements_list <-lapply(seq_along(gene_lists$gene_list), function(i) {# Current target list target_list <- list_of_lists[[i]]# All other lists except the current one other_lists <- list_of_lists[-i]# Find unique elements for the current listfind_unique_elements(target_list, other_lists)})# Assign names to the resulting list of listsnames(unique_elements_list) <-names(list_of_lists)
'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns
building nested credible set table
Joining with `by = join_by(geneSetIdx)`
Joining with `by = join_by(geneSetIdx)`
Joining with `by = join_by(alpha, component, geneSet)`
Joining with `by = join_by(geneSet)`
Joining with `by = join_by(gene)`
Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
# Was:
data %>% select(component_columns)
# Now:
data %>% select(all_of(component_columns))
See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
# Was:
data %>% select(gene_set_columns)
# Now:
data %>% select(all_of(gene_set_columns))
See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
Warning: Unknown or uninitialised column: `details`.
Unknown or uninitialised column: `details`.
Unknown or uninitialised column: `details`.